# region 针对id为24的门，进行了关门的最大位移统计分析
# 使用逻辑回归模型训练该门体的健康评估模型
# 训练集为2021年4月当中的13天(6,8,10,13,14,15,16,18,19,20,26,27,29)的数据
# 测试集为20210501,20210502,20210430这3天的数据
# 训练集和测试集，健康和亚健康，的数量比，控制在了大约4：1的样子
# print(len(np.where(y_train==0)[0]),len(np.where(y_train==1)[0]),count_train)    # 12 2009 2021
# print(len(np.where(y_test==0)[0]),len(np.where(y_test==1)[0]),count_test)       # 3  506  509
# 设定五个位移位置：200，400，600，800，smax
from __future__ import division
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymysql
from numpy.ma.core import count
from scipy import integrate
from scipy.integrate._ivp.radau import P
from scipy.optimize import curve_fit
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import precision_score, recall_score
from sklearn.model_selection import GridSearchCV
warnings.filterwarnings("ignore")
from matplotlib.colors import ListedColormap

# endregion

# region 函数定义
def getv(df):
    """
    对每一条dataframe进行处理并存储到一维数组v
    :param df: 从数据库读取到的数据帧
    :return: 速度数组v
    """
    try:
        # 定义data数组复制一条datas
        data = np.arange(220)  
        for i in range(220):
            data[i] = df[0][i]
        # speed数组的长度是每条datas的一半
        v = np.arange(110)  
        a = 0
        for i in range(110):
            v[i] = data[a]*256+data[a+1]     # 转换公式 result = 256 * x + y x：第一个字节，Y_train:第二个字节
            a = a + 2
        return v
    except Exception as e:
        print(type(e),e)


def gets(count,v):
    """
    通过速度求位移
    :param count: 读取到的datas的数量
    :param v: 处理datas后得到的速度
    :return: 位移数组s
    """
    try:
        # 定义位移数组
        s = np.zeros((count,110))
        for i in range(1,count):      # 中间值速度乘以采样时间间隔得到位移
            if(i<110):
                s[0][i] = s[0][i-1] + (v[0][i-1]+v[0][i]) / 2.0 * 0.045
            for j in range(1,110):
                s[i][j] = s[i][j-1] + (v[i][j-1]+v[i][j]) / 2 * 0.045 
        return s
    except Exception as e:
        print(type(e),e)


def getsmax(count,s):
    """
    通过位移数组取最大位移
    :param count: 读取到的datas的数量
    :param s: 得到的位移数组
    :return: 最大位移数组smax
    """
    try:
        smax = np.zeros(count)
        for i in range(count):
            smax[i] = s[i][len(s[i])-1]
        return smax   
    except Exception as e:
        print(type(e),e)


def gety(count,smax,mean,std):
    """
    通过最大位移的均值和方差判断特征y为健康--1或亚健康--0
    :param count: 读取到的datas的数量
    :param smax: 最大位移数组
    :param mean: 最大位移的均值
    :param std: 最大位移的方差
    :return: 特征y, 健康--1, 亚健康--0
    """
    try:
        y = np.ones(count)
        for i in range(count):
            if smax[i] > mean+std*3 or smax[i] < mean-std_train*3:
                y[i] = 0
        return y
    except Exception as e:
        print(type(e),e)


def getsmaxtime(s,smax):
    """
    求解首次到达最大位移所花费的时间
    :param s: 此次的位移数组
    :param smax: 此次的最大位移
    :return: 首次到达最大位移所花费的时间
    """
    try:
        t = 0.0
        n = len(s)
        i = 0
        while i < n:
            if s[i] >= smax:
                t = i*0.045
                break
            i += 1
        return t
    except Exception as e:
        print(type(e),e)

def getx(count,s):
    """
    通过位移求出输入变量x，x为门体从指定位置到达指定位置所花费的时间
    设定五个位移位置：200，400，600，800，smax
    :param count: 读取到的datas的数量
    :param s: 位移数组
    :return: 输入变量数组x
    """
    try:
        x = np.zeros((count,5))
        for i in range(count):
            if s[i][109]<200:
                continue
            elif s[i][109]>=200 and s[i][109]<400:
                x[i][0] = np.interp(200,s[i],t)
            elif s[i][109]>=400 and s[i][109]<600:
                x[i][0] = np.interp(200,s[i],t)
                x[i][1] = np.interp(400,s[i],t) - x[i][0]
            elif s[i][109]>=600 and s[i][109]<800:
                x[i][0] = np.interp(200,s[i],t)
                x[i][1] = np.interp(400,s[i],t) - x[i][0]
                x[i][2] = np.interp(600,s[i],t) - x[i][0] - x[i][1]
            elif s[i][109]==800:
                x[i][0] = np.interp(200,s[i],t)
                x[i][1] = np.interp(400,s[i],t) - x[i][0]
                x[i][2] = np.interp(600,s[i],t) - x[i][0] - x[i][1]
                x[i][3] = np.interp(800,s[i],t) - x[i][0] - x[i][1] - x[i][2]
            elif s[i][109]>800:
                x[i][0] = np.interp(200,s[i],t)
                x[i][1] = np.interp(400,s[i],t) - x[i][0]
                x[i][2] = np.interp(600,s[i],t) - x[i][0] - x[i][1]
                x[i][3] = np.interp(800,s[i],t) - x[i][0] - x[i][1] - x[i][2]
                x[i][4] = getsmaxtime(s[i],s[i][109]) - x[i][0] - x[i][1] - x[i][2] - x[i][3]
        return x
    except Exception as e:
        print(type(e),e)
def getblr(lr,x,y):
    """
    训练获取最优参数，获取调参后的逻辑回归模型blr
    :param lr: 未经过调参的逻辑回归函数
    :param x: 训练集的x
    :param y: 训练集的y
    :return: 调参后的逻辑回归模型blr和最优参数组合best_params
    """
    try:
        params = {'C':[0.0001, 1, 100, 1000],
            'max_iter':[1, 10, 100, 500],
            'class_weight':['balanced', None],
            'solver':['liblinear','sag','lbfgs','newton-cg']
            }
        clf = GridSearchCV(lr, param_grid=params, cv=10)
        clf.fit(x,y)
        best_params = clf.best_params_  
        blr = LogisticRegression(**best_params)
        blr.fit(x,y)
        return blr,best_params
    except Exception as e:
        print(type(e),e)
def dtg(s,ts):
    """
    绘制位移时间图
    :param s: 位移数组
    :param ts: 到达指定位移所花费的时间
    :return: none
    """
    try:
        plt.plot(t,s)
        plt.axhline(200, color='r', linestyle='--', label='200')
        plt.axhline(400, color='r', linestyle='--', label='400')
        plt.axhline(600, color='r', linestyle='--', label='600')
        plt.axhline(800, color='r', linestyle='--', label='800')
        plt.axhline(s[len(s)-1], color='g', linestyle='--', label='max')
        plt.axvline(ts[0], color='b', linestyle='--', label='t1')
        plt.axvline(ts[0]+ts[1], color='b', linestyle='--', label='t2')
        plt.axvline(ts[0]+ts[1]+ts[2], color='b', linestyle='--', label='t3')
        plt.axvline(ts[0]+ts[1]+ts[2]+ts[3], color='b', linestyle='--', label='t4')
        plt.axvline(ts[0]+ts[1]+ts[2]+ts[3]+ts[4], color='b', linestyle='--', label='t5')
        plt.legend()
        plt.xlabel('time')
        plt.ylabel('displacement')
        plt.title('Displacement-time graph')
        plt.show()
    except Exception as e:
        print(type(e),e)
def mdtg(count,smax,mean,std):
    """
    绘制最大位移时间图
    :param count: 最大位移数组长度(最大位移数量)
    :param smax: 最大位移数组
    :param mean: 最大位移数组的均值
    :param std: 最大位移数组的标准差
    :return: none
    """
    try:
        x = np.arange(count)
        plt.plot(x,smax,marker='.')
        plt.axhline(mean, color='r', linestyle='--', label='mean')
        plt.axhline(mean+std, color='r', linestyle='--', label='mean+std')
        plt.axhline(mean-std, color='r', linestyle='--', label='mean-std')
        plt.axhline(mean+std*2, color='r', linestyle='--', label='mean+std*2')
        plt.axhline(mean-std*2, color='r', linestyle='--', label='mean-std*2')
        plt.axhline(mean+std*3, color='r', linestyle='--', label='mean+std*3')
        plt.axhline(mean-std*3, color='r', linestyle='--', label='mean-std*3')
        plt.legend(loc='lower left')
        plt.xlabel('time')
        plt.ylabel('maxdisplacement')
        plt.title('Maxdisplacement-time Graph')
        plt.show()
    except Exception as e:
        print(type(e),e)
def slrg():
    """
    用最大位移和刚达到最大位移的时间作为特征，绘画逻辑回归分类图
    :return: none
    """
    try:
        smaxtime_train,smaxtime_test = np.zeros(count_train),np.zeros(count_test)
        for i in range(count_train):
            smaxtime_train[i] = getsmaxtime(s_train[i],smax_train[i])
        for i in range(count_test):
            smaxtime_test[i] = getsmaxtime(s_test[i],smax_test[i])
        X_train,X_test = np.vstack((smaxtime_train,smax_train)).T , np.vstack((smaxtime_test,smax_test)).T
        slr = LogisticRegression()
        slr.fit(X_train,y_train)
        X_set,y_set = X_test,y_test
        X1,X2=np.meshgrid(np.arange(start=X_set[:,0].min()-1, stop=X_set[:, 0].max()+1, step=0.01),
                        np.arange(start=X_set[:,1].min()-1, stop=X_set[:,1].max()+1, step=0.01))
        plt.contourf(X1, X2, slr.predict(np.array([X1.ravel(),X2.ravel()]).T).reshape(X1.shape),
                    alpha = 0.75, cmap = ListedColormap(('#C65749', '#338DFF')) )
        plt.xlim(X1.min(),X1.max())
        plt.ylim(X2.min(),X2.max())
        for i,j in enumerate(np.unique(y_set)):
            plt.scatter(X_set[y_set==j,0],X_set[y_set==j,1],
                        color = ListedColormap(('red', 'blue'))(i), label=j)
        plt.title(' LOGISTIC(Test set)')
        plt.xlabel(' smaxtime')
        plt.ylabel(' smax')
        plt.legend()
        plt.show()
    except Exception as e:
        print(type(e),e)

# endregion

# region 获取数据，处理得到训练集和测试集
# 连接数据库，获取数据
db = pymysql.connect(host='10.10.62.36', user='root', password='root',database='stationdb', charset='utf8')
sql_train = " SELECT datas FROM db_closedoorspeedcurve_20210406 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210408 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210410 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210413 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210414 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210415 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210416 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210418 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210419 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210420 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210426 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210427 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210429 WHERE Id='24' "
sql_test = "SELECT datas FROM db_closedoorspeedcurve_20210501 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210502 WHERE Id='24' UNION ALL \
        SELECT datas FROM db_closedoorspeedcurve_20210430 WHERE Id='24' "
df_train,df_test = pd.read_sql_query(sql_train,db),pd.read_sql_query(sql_test,db)
# 记录从数据库得到的datas的数量
count_train,count_test = len(df_train['datas']),len(df_test['datas'])
# 对每个datas处理得到想要的速度数据并存储到二维数组
v_train,v_test = df_train.apply(getv,axis=1),df_test.apply(getv,axis=1)
# 时间从0.0开始，每0.045s采样一次，共采样了110次，到4.905s结束
t = np.arange(0.0,4.905,0.045)
# 定义二维数组存放位移
s_train,s_test = gets(count_train,v_train),gets(count_test,v_test)
# s_train,s_test = gets2(count_train,v_train),gets2(count_test,v_test)
# 记录每条数据的最大位移
smax_train,smax_test = getsmax(count_train,s_train),getsmax(count_test,s_test)
# 计算均值和标准差
mean_train,mean_test,std_train,std_test = np.mean(smax_train),np.mean(smax_test),np.std(smax_train),np.std(smax_test)
# 记录是否健康这一特征--健康为1，亚健康为0，对于超出3std的数据判定为亚健康
y_train,y_test = gety(count_train,smax_train,mean_train,std_train),gety(count_test,smax_test,mean_test,std_test)
# 记录每条数据的五个特定时间作为输入变量，达不到指定位移的话将该时间段设为0
x_train,x_test = getx(count_train,s_train),getx(count_test,s_test)
# 断开和数据库的连接
db.close()
# endregion

# region 用sklearn进行逻辑回归
lr = LogisticRegression()
lr.fit(x_train,y_train)
lrpred = lr.predict(x_test)
blr,best_params = getblr(lr,x_train,y_train)
blrpred = blr.predict(x_test)
#  print("未经过调参的模型的精确率和召回率：",precision_score(y_test,lrpred),recall_score(y_test,lrpred))
# print("训练得出的最优参数组合：", best_params)
# print("经过调参优化的模型的精确率和召回率：",precision_score(y_test,blrpred),recall_score(y_test,blrpred))
# dtg(s_train[0],x_train[0])
dtg(s_test[243],x_test[243])
dtg(s_test[422],x_test[422])
dtg(s_test[469],x_test[469])
# print(np.where(y_test==0)[0])   # [243,422,469]
def vtg(t,v):
    plt.plot(t,v)
    plt.show()
vtg(t,v_test[243])
vtg(t,v_test[422])
vtg(t,v_test[469])
# mdtg(count_test,smax_test,mean_test,std_test)
# slrg()
# endregion

# region 对由速度求位移的三种方法进行研究
# # 通过数据曲线拟合得到一个多项式,用该多项式求积分
# func = np.poly1d(np.polyfit(t, v_train[0], 11)) # 936.7213886571899
# func = np.poly1d(np.polyfit(t, v_train[0], 12)) # 936.4719344810686
# def curvefitgraph(func):
#     """
#     画图查看曲线的拟合程度
#     :return: none
#     """
#     try:
#         yvals = func(t) 
#         plt.plot(t, v_train[0], label='original values')
#         plt.plot(t, yvals, label='polyfit values')
#         plt.title('curve_fit')
#         plt.show()
#     except Exception as e:
#       print(type(e),e)
# def gets2(count,v):
#     """
#     对速度积分求位移
#     :param count: 读取到的datas的数量
#     :param v: 处理datas后得到的速度
#     :return: 位移数组s
#     """
#     try:
#         # 定义位移数组
#         s = np.zeros((count,110))
#         for i in range(count):
#             func = np.poly1d(np.polyfit(t, v[i], 11)) 
#             for j in range(110):
#                 s[i][j], err = integrate.quad(func, 0.0, 0.045*(j+1))
#         return s
#     except Exception as e:
#         print(type(e),e)
# curvefitgraph(func) # 查看曲线拟合程度
# # 通过更改积分上限可以求出不同时刻的位移
# s1, err = integrate.quad(func, 0.0, 4.905)
# print("用曲线拟合得到的十一次多项式求的积分：",s1)
# # 在无法确认积分函数的情况下，给出一些序列也可做积分。但是这种方法不能求出每个时刻的位移，只能一次求出最大位移
# s2 = integrate.trapz(v_train[0], t)
# print("直接用训练集的序列求的积分：",s2)
# # 与中值法进行对比
# print("之前用的中值法求的位移：",smax_train[0])

# endregion

